The two-loop self-energy for the ground state of medium- Z hydrogen-like ions 
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The two-loop self-energy correction to the ground state Lamb shift is calculated for hydrogen-like ions with 
the nuclear charge Z — 10 ... 30 without any expansion in the binding field of the nucleus. A calculational 
technique is reported for treatment of Feynman diagrams in the mixed coordinate-momentum representation, 
which yields significant improvement in numerical accuracy as compared to previous results. An extrapolation 
of the all-order numerical data towards Z — 1 yields a result for the higher-order remainder function for hydro- 
gen. The previously reported disagreement between the all-order and the perturbative approaches is reduced to 
the marginal agreement. 
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Calculations of the two-loop QED corrections to all orders 
in the binding-strength parameter of the nucleus Za (Z is 
the nuclear charge and a is the fine structure constant) were 
motivated almost two decades ago by progress injhe experi- 
mental spectroscopy of heavy lithium-like ions |1]. To date, 
experimental investigations of the 2pj-2s transitions in high- 
and medium- Z lithium-Uke ions are sensitive to the two-loop 
QED effects on the level of about 10% f^, All-order 
calculations of the two-loop Lamb shift were recently ac- 
complished in Ref. |4] for high-Z hydrogen-like ions and in 
Ref. [5^ for high-Z lithium-like ions. 

The ratio of the QED effects to the binding energy scales 
as Z^, so that the relative QED contribution gets smaller for 
lighter ions. However, the experimental precision is also bet- 
ter there and the two-loop QED effects in light systems have 
long being observed. The best studied case is atomic hydro- 
gen, whose spectroscopy is nowadays carried out with an ac- 
curacy of a few parts in 10^"* Q]. 

Until recently, calculations of the QED effects in hydrogen 
relied on the approach perturbative in the binding-strength pa- 
rameter Za. Technical difficulties of this approach, however, 
grow rapidly with the increase of the order of the perturba- 
tive expansion. The state of the art of such calculations is the 
evaluation of the dominant part of the (Za)^ correction 
10 Isl Ist]- An alternative way is to perform a numerical cal- 
culation to all orders in Za and to identify the higher-order 
remainder by subtracting the known low-order terms from the 
all-order results. The problems on this way are, first, signif- 
icant internal cancellations in numerical calculations, which 
grow as Z decreases, and second, additional losses of ac- 
curacy occuring when the higher-order remainder is inferred 
from the all-order results. 

The first attempt at the evaluation of the all-order remain- 
der for the ground state of hydrogen was made in Ref. ifToll . 
In that work, a numerical calculation of the two-loop self- 
energy correction was reported for Z > 10. This correction 
is expected to give the dominant contribution to the two-loop 
remainder An extrapolation towards Z = 1 performed in 
Ref. ifToll yielded a result for the higher-order remainder ap- 
proximately twice as large as the estimate based on the an- 
alytical calculations |8]. This disagreement is presently the 
main source of the theoretical uncertainty of the ground-state 
Lamb shift in hydrogen and influences the values of the Ryd- 



berg constant and the proton charge radius obtained from the 
hydrogen spectroscopic data lITlll . 

This investigation presents an attempt to resolve the dis- 
agreement between the numerical and analytical approaches 
by improving the calculational accuracy of the all-order re- 
sults. To this end, we develop a scheme for the evaluation of 
Feynman diagrams in the mixed coordinate-momentum repre- 
sentation (the corresponding part of the two-loop self-energy 
is conventionally termed as the P term). This is the most non- 
trivial part of the evaluation of the two-loop self-energy cor- 
rection, as it has no analog in the one-loop calculations. 

For the first time the P term was calculated in Ref. lfl2ll 
with help of a finite basis set representation of the spectrum of 
the Dirac equation. Later investigations |4, 5, proceeded 
along the same way with adopting an improved (dual kinet- 
ically ballanced) basis set lUSll . The main problem of the 
basis-set approach is a relatively slow convergence with re- 
spect to the number of basis functions. In order to overcome 
this limitation, in this work we employ the analytical repre- 
sentation of the Dirac-Coulomb Green function (DCGF) in 
terms of the Whittaker functions. As a result, the numerical 
accuracy of the P term is improved by more than an order of 
magnitude, the error now being mainly due to the termination 
of the partial-wave expansion. 

In the present investigation, we perform a reevaluation of 
the two-loop self-energy correction for the ground state of 
hydrogen-like ions with the nuclear charge numbers in the in- 
terval Z ^ 10 ... 30. The calculation of the P term is carried 
out with the technique developed in this work. The other parts 
of the correction are evaluated by the methods described pre- 
viously [14] but with the increased number of partial waves 
included and with denser integration grids. The nonpertur- 
bative remainder incorporating terms of order (Za)^ and 
higher is inferred from the numerical results and extrapolated 
toward Z — 1. 

The two-loop self-energy correction is shown in Fig.[T] In 
this work, we concentrate on the part of the correction that 
is treated in the mixed coordinate-momentum representation 
and is referred to as the P term. The corresponding Feynman 
diagrams are shown in Fig. |2] They arise from the diagrams 
in Fig. [T] when the bound-electron propagators are expanded 
in terms of the interaction with the binding field. The dis- 
tinct feature of the diagrams contributing to the P term is that 
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FIG. 1 : The two-loop self-energy correction. The double line rep- 
resents the electron propagator in the Coulomb binding field of the 
nucleus. 
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FIG. 2: The P term. The single line represents the free electron 
propagator. The dashed line with a cross indicates the interaction 
with the Coulomb field of the nucleus. 



ultraviolet divergences in them originate from the one-loop 
subgraphs only. These subgraphs are covariantly regularized 
and calculated in the momentum space, whereas the remain- 
ing part of the diagrams does not need any regularization and 
is treated in the coordinate space. 

In order to illustrate the calculational technique used for 
the evaluation of the P term, we consider one of the simplest 
diagrams in Fig. |2] graph (b). Its contribution can be written 
as 
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(1) 



where E = Ea—oj, Sa is the reference-state energy, D{u!, X14) 
is the scalar part of the photon propagator in the Feynman 
gauge defined by X14) = .g'^" D{lu, xu), is the 

renormalized free self-energy operator (see, e.g., Ref. lITill 
for its definition and evaluation), Gv{£,xi,p) is the Fourier 
transform of the product of the DCGF G(e, xi,X2) and the 
Coulomb potential Vc, 



Gv{e,xi,p)^ J dx2 e'P=^^ G{e,Xi,X2)Vc{x2) 



(2) 



the function Gy is given by an analogous expression with 
the DCGF substituted by the free Dirac Green function G^°\ 
and G'"' (E, p, X4^) is the Fourier transform of the free Dirac 
Green function. 

The main problem of the evaluation of Eq. ([T]! is the absence 
of a convenient representation for the function Gy(e, Xi,p). 



Because of this, we have to start with the the coordinate-space 
representation of the Green function and perform the Fourier 
transfomation numerically. The DCGF is known in terms of 
the partial-wave expansion over the relativistic angular param- 
eter K, with the radial part given by [15,1 



GK{e,xi,X2) = 



-0°(£,a:i)(/)^''(£,X2), forxi<a;2, 
-(j)'^{e,xi)(j)f {e,x2) , forx2<xi, 

(3) 

where (j)^ and (j)"^ are the two-component solutions of the ra- 
dial Dirac equation regular at the origin and at the infinity, 
respectively, and normalized in such a way that their Wron- 
skian is unity. The Fourier transform of the DCGF over the 
second radial variable can be represented as 



Gk(£,xi,p) 



(4) 



with 



V'°(£,p;xi) = 47r 
and 

ij^{e,p;xi) ^4tt 



dx2 X2 



dX2 x\ 



3l{pX2)4>l^+{£-,X2) 

^ir(pa;2)0°,_(£,a;2) 



ji{px2)(l)'^^+{e,X2) \ 

J^jjipX2)<P^^^{£,X2) J 

(6) 



where (/)° _|_ and 



denote the upper and lower components 
respectively, / = |k + 1/2| — 1/2 and I = 



of (jr^ and c 
\k- 1/2| - 1/2. 

The integration over X2 in the functions ilji^{£,p; xi) has to 
be performed numerically. The problems here are that (i) the 
integration interval depends on xi and (ii) the integrand con- 
tains the spherical Bessel function which oscillates rapidly in 
the high-momenta region. Clearly, a straightforward use of 
Eqs. ^ and (|6]l would lead to a re-evaluation of the integral 
for each new value of xi, making the calculation prohibitively 
expensive. One can observe, however, that if the function 
4>k{s,p; x) is known for a particular set of e, p, and x, then 
the evaluation of ■ipi^{e,p; x') can be done by computing the 
Bessel transform integral over the interval {x,x') only. So, 
introducing an ordered radial grid {xi}, one can store the set 
of values {V'k(£,p; Xi)} by performing just one Bessel trans- 
form over the interval (0, 00). This shows that for a fixed 
values of e and p, the integration over xi can be performed 
without a recalculation of the Bessel transform integral. Still, 
the evaluation of the functions ipk was one of the most prob- 
lematic parts of the computation since a controllable accuracy 
was required for momenta as high as p = 10^. 

The next problem to be solved in the numerical evaluation 
of Eq. ([T]l is that the free Dirac Green function G*-*^-* {E, p, x^) 
contains a spherical Bessel function jhipxA) and thus is 
highly oscillating too in the high-momenta region. It can 
be observed that the radial integration over xa resembles the 
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FIG. 3: The higher-order remainder function of the two-loop self- 
energy correction. The cross on the y-axis indicates the analytical 
result. 



Bessel transform integral over X2 and so can be efficiently cal- 
culated by introducing analogs of the functions tjjk- Details of 
the numerical procedure will be published elsewhere. 

The results of our numerical evaluation of the two-loop self- 
energy correction to the ground-state Lamb shift of middle-Z 
hydrogen-like ions are presented in Table |I] They are con- 
sistent with but improve upon the data obtained previously 

MM- 

In the present investigation we are concerned with the 
higher-order remainder function that incorporates contribu- 
tions of all orders starting with a'^{Za)^ and is denoted as 
^SESE- obtained from the two-loop self-energy correc- 
tion Ai?sESE by separating out the first terms of its Za ex- 
pansion, 



>50 



A£;sESE = m |s4o + (^a)B, 



(7) 



where L = \n[{Za)^'^] and the expansion of the remainder 
starts with a constant, Ggggg(2') = Bqq + Za(..A. The 
results for the expansion coefficients (see Refs. Q H HI [Ull 
and references therein) are: B40 — 1.409244, B50 = 
-24.2668(31), Bgs = -8/27, ^62 = 16/ 27 - (16/9) ln2, 
Bei = 48.388913, and Bqq = -61.6(9.2). The remain- 
der function infeiTed from our numerical results is plotted in 
Fig.S 

In order to obtain a value of the remainder function for hy- 
drogen, we have to extrapolate the numerical data obtained 
for Z > 10 towards Z = 1. For this we use a variant of the 



procedure first employed in Ref . (vf] . The extrapolated value 
of Ggggg(Z = 1) is obtained in two steps. Fkst, we apply an 
(exact) linear fit to each pair of two consecutive points from 
our data set and store the resulting values at Z = 1. Sec- 
ond, we perform a global parabolic least-squares fit to the set 
of data obtained on the first step and take the fitted value at 
Z = 1 as a final result. Similar procedure applied to the de- 
termination of the B^o coefficient reproduces the known ana- 
lytical result with the accuracy of about 1%. For comparison, 
a global polynomial fit yields a result for the B^q coefficient 
accurate within 5% only. 

When applied to the remainder function Gg^g^{Z), the ex- 
trapolation procedure described above gives 



G§ese(^ = 1) = -86(15) 



(8) 



This value is significantly higher than what would be expected 
simply from a polynomial fitting (of about —105), which sug- 
gests the presence of unusually large logarithmic contribu- 
tions to order oHZaY (logarithms of up to the second power 
are anticipated H). The value (O is also higher than but con- 
sistent with the 2005 result of -127 (42) |10]. The shift of 
the central value is due to two reasons. First, the analytical 
value of the Sgi coefficient was changed during this time by 
SBqi = —1.4494 . . . |9], thus pushing the higher-remainder 
higher up. Second, the improved numerical accuracy of the 
present calculation and the increased number of data points 
allowed us to identify the upward trend in the numerical data, 
which influenced the extrapolated value considerably. The 
present result ([8]) is marginally consistent with the analytical 
value of Seo = -62(9). 

Finally, we account for the contribution from the diagrams 
with closed fermion loops calculated recently in Ref. 1. 1 81I . 
The total value of the remainder function for hydrogen is 



G^ ° {Z = 1) = -86 (15) - 15 (2) = -101 (15) 



(9) 



To conclude, the present investigation reports a tech- 
nique for the evaluation of Feynman diagrams in the mixed 
coordinate-momentum representation, which allows one to 
significantly improve the numerical accuracy. A complete 
recalculation of the diagrams of the two-loop self-energy is 
presented for the ground state of hydrogen-like ions with the 
nuclear charge number Z — 10 — 30. The higher-order (in 
Za) remainder function is inferred from numerical all-order 
results and extrapolated towards Z = 1. The extrapolated 
value of the higher-order remainder function is in marginal 
agreement with the analytical result obtained within the per- 
turbative approach. 

The work reported in this paper was supported by the "Dy- 
nasty" foundation. 
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